rm(list=ls())
library(ggplot2)
library(readstata13)
library(rgdal)
library(scales)
library(ggpubr)

options(scipen=999)
setwd("H:/Stata/Data/temp/")

data <- read.dta13("./toplot.dta") 



###################################################
#Figure A1: Variation in Treatment Proportion
###################################################

ggplot(data, aes(treat_change_pop,treat_total_2C_pop, color=treat_change_pop)) +
  geom_point(shape=19, size=5.5, alpha=.7) + 
  theme_minimal() + 
    scale_color_gradient2(low="darkblue", mid="skyblue" , high="aquamarine") + 
    theme(legend.position = "none") + 
  xlab("Change in Treatment") + ylab("Mean % Treated") + 
  scale_y_continuous(breaks=c(0,0.005, 0.010)) 
ggsave("../../Output/toexport_figures/figa1.pdf")

###################################################
#Figure 2: Left Vote Share and Change in \% Treated
###################################################

ggplot(data, aes(treat_change_pop,left_pct, color=treat_change_pop)) +
  geom_point(shape=19, size=5.5, alpha=.7) + 
  theme_minimal() + 
  scale_color_gradient2(low="darkblue", mid="skyblue" , high="aquamarine") + 
  theme(legend.position = "none") + 
  xlab("Change in Treatment") + ylab("Mean Left Support") +
  geom_smooth(method="lm", se=FALSE, color="red", fill="grey100", lwd=1.2)
ggsave("../../Output/toexport_figures/fig2a.pdf")

ggplot(data, aes(treat_change_pop,left_change, color=treat_change_pop)) +
  geom_point(shape=19, size=5.5, alpha=.7) + 
  theme_minimal() + 
  scale_color_gradient2(low="darkblue", mid="skyblue" , high="aquamarine") + 
  theme(legend.position = "none") + 
  xlab("Change in Treatment") + ylab("Change in Left Support") +
  geom_smooth(method="lm", se=FALSE, color="red", fill="grey100", lwd=1.2)
ggsave("../../Output/toexport_figures/fig2b.pdf")


###################################################
#Figure A2: Left Vote Share and $\Delta \%$ Treated Within West and East
###################################################

ggplot(data[data$west_provincie==1,], aes(treat_change_pop,left_change, color=treat_change_pop)) +
  geom_point(shape=19, size=5.5, alpha=.7) + 
  theme_minimal() + 
  scale_color_gradient2(low="darkblue", mid="skyblue" , high="aquamarine") + 
  theme(legend.position = "none") + 
  xlab("Change in Treatment") + ylab("Change in Left Support") +
  geom_smooth(method="lm", se=FALSE, color="red", fill="grey100", lwd=1.2) + 
  ggtitle("West Municipalities")
ggsave("../../Output/toexport_figures/figa2a.pdf")

ggplot(data[data$west_provincie==0,], aes(treat_change_pop,left_change, color=treat_change_pop)) +
  geom_point(shape=19, size=5.5, alpha=.7) + 
  theme_minimal() + 
  scale_color_gradient2(low="darkblue", mid="skyblue" , high="aquamarine") + 
  theme(legend.position = "none") + 
  xlab("Change in Treatment") + ylab("Change in Left Support") +
  geom_smooth(method="lm", se=FALSE, color="red", fill="grey100", lwd=1.2) + 
  ggtitle("East Municipalities")

ggsave("../../Output/toexport_figures/figa2b.pdf")

###################################################
#Figure A3: Average \textbf{Levels} and \textbf{Change} of Treated
###################################################

myshp <- "H:/Stata/Municipalities and provinces/Shapefiles/Gemeenten_2017__CBS_Wijk_en_Buurtkaart.shp"
myfile <- readOGR(myshp)


df <- fortify(myfile)
df$id <- as.numeric(as.character(df$id))
df$id <- df$id + 1 ## to match below, ids from 1 to 388 
summary(df$id)



munid <- data.frame(myfile$gm_code, myfile$objectid)
names(munid) <- c("munname", "id")
munid$id <- as.numeric(as.character(munid$id))
summary(munid$id)

merge <- merge(df, munid, by="id")

muntreat <- read.dta13("H:/Stata/Do-files/municipality/maps/mapsmun_map_data.dta")

data.final <- merge(merge, muntreat, by="munname")

vars <- c("munname", "id", "long", "lat", "munid", 
          "meanMUN_treat", "meanMUN_left", "group", 
          "meanMUN_treatchange", "meanMUN_leftchange", 
          "meanMUN_treataggch", "meanMUN_treatagg")
data.final <- data.final[vars]


ggplot(data.final, aes(x=long, y=lat, group=group)) + geom_path() + 
  coord_map("mercator") + geom_polygon(aes(fill=(meanMUN_treatagg)), color="white") +
  theme_minimal() + 
  xlab(" ") + ylab(" ") +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text = element_blank(), 
        legend.key.size = unit(.3,"cm"),
        legend.text = element_text(size=6),
        legend.title = element_text(size=6)) +
  scale_fill_gradientn(colors = c("skyblue", "blue", "darkblue"), 
                       values=scales::rescale(c(0, 250, 500, 1000)), 
                       name="Treated Level") 
ggsave(filename = "./toexport_figures/figa3a.pdf")


ggplot(data.final, aes(x=long, y=lat, group=group)) + geom_path() + 
  coord_map("mercator") + geom_polygon(aes(fill=(meanMUN_treataggch)), color="white") +
  theme_minimal() + 
  xlab(" ") + ylab(" ") +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text = element_blank(),
        legend.key.size = unit(.3,"cm"),
        legend.text = element_text(size=6),
        legend.title = element_text(size=6)) + 
  scale_fill_gradientn(colors = c("tomato", "gray70", "skyblue"), 
                       values=scales::rescale(c(-50, -20, 
                                                0, 10, 20)), 
                       name="Treated Change")
ggsave(filename = "./toexport_figures/figa3b.pdf")

###################################################
#Figure 1: Regions Affected by the Dutch Winter Hunger
###################################################

myshp <- "NLD_adm1.shp"
myfile<-readOGR(myshp)

df <- fortify(myfile)
df$id <- as.numeric(df$id)
summary(df$id) ## 14 regions 

names <- read.delim(file="NLD_adm1.csv", sep=",")
head(names)
names$ID_1  <- names$ID_1 -1 ## this can be merged to df to identify ID 

#treated provinces: North and South Holland and Utrecht.
df$famine <- 0 
df$famine <- ifelse(df$id==8, 1, df$famine) ## north holland
df$famine <- ifelse(df$id==13, 1, df$famine) ## south holland
df$famine <- ifelse(df$id==10, 1, df$famine) ## utrecth

df$famine <- as.factor(df$famine)

gg <- ggplot(df, aes(x=long, y=lat, group=group)) + geom_path() + coord_map("mercator") + 
  geom_polygon(aes(fill=famine),  color = "white") + 
  theme_minimal() + 
  xlab(" ") + ylab(" ") + 
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text = element_blank()) + 
  scale_fill_manual(values=c("grey60", "goldenrod1"),  breaks=c("0", "1"),
                    labels=c("Control", "Treated"), name="Famine Status") 

gg

ggsave(gg, file="fig1.pdf")


###################################################
#Figure 3: Left Support in West Pre Famine
###################################################

# Socio Effects
data2 <- read.dta13("./regsave1.dta")
data2$varorder <- c(1,2,3,4)
ggplot(data2, aes(as.factor(varorder), coef)) + 
  geom_pointrange(aes(ymin=coef-1.96*stderr, ymax=coef+1.96*stderr), shape=16, size=1.2, 
                  color="skyblue") + 
  geom_hline(yintercept = 0, lwd=1.2, color="red") +
  theme_minimal() + ylab("Effect on Left Vote Share") + xlab(" ") +  
  theme(axis.title = element_text(size=12),
        axis.text = element_text(size=12)) + 
  scale_y_continuous(breaks=c(-.02, -0.01, 0, 0.01, 0.02), limits=c(-.02, 0.02)) + 
  scale_x_discrete(labels=c("Pre W 10yo", "Pre W 5yo", "Pre W 2yo", "Treat"))
ggsave("../../Output/toexport_figures/fig3.pdf")

###################################################
#Figure 4: Left Support in West Post Famine
###################################################

data3 <- read.dta13("./regsave2.dta")
data3$varorder <- c(1,2,3,4)
ggplot(data3, aes(as.factor(varorder), coef)) + 
  geom_pointrange(aes(ymin=coef-1.96*stderr, ymax=coef+1.96*stderr), shape=16, size=1.2, 
                  color="skyblue") + 
  geom_hline(yintercept = 0, lwd=1.2, color="red") +
  theme_minimal() + ylab("Effect on Left Vote Share") + xlab(" ") +  
  theme(axis.title = element_text(size=12),
        axis.text = element_text(size=12)) + 
  scale_y_continuous(breaks=c(-.02, -0.01, 0, 0.01, 0.02), limits=c(-.02, 0.02)) + 
  scale_x_discrete(labels=c("Treat", "Post W 2yo", "Post W 5yo", "Post W 10yo"))
ggsave("../../Output/toexport_figures/fig4.pdf")

###################################################
#Figure 5: Treated Movers and Municipality Characteristics
###################################################

dt <- read.dta13("./regsave_sorting_ncontrols.dta")
dt$varorder <- c(5,4,3,2,1)
dt$varorder2 <- c(1,2,3,4,5)

ggplot(dt, aes(as.factor(varorder2), coef)) + 
  geom_pointrange(aes(ymin=coef-1.96*stderr, ymax=coef+1.96*stderr), 
                  shape=16, size=1.2,  color="skyblue") + 
  geom_hline(yintercept = 0, lwd=1.2, color="red") +
  theme_minimal() + ylab("Effect of %Treated Mover On...") + xlab(" ") +  
  theme(axis.title = element_text(size=12),
        axis.text = element_text(size=12)) + 
  scale_y_continuous(breaks=c(-.05, 0, 0.05), limits=c(-.05, 0.05)) + 
  scale_x_discrete(labels=c("Log Social Spending","Log Education Spending",
                            "Log Health Spending", "Female Pop", "Log Home Price"))+
  coord_flip()
ggsave("../../Output/toexport_figures/fig5.pdf")
